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ABSTRACT 


Specification of the mathematical equations required to define the dynamic response of a linear con- 
tinuous plant, subject to sampled data control, is complicated by the fact that the digital components 
of the control system cannot be modeled via linear ordinary differential equations. This complication 
can be overcome by introducing two new mathematical operations; namely, the operation of zero 
order hold and digital delay. It is shown herein that by direct utilization of these operations, a set of 
linear mixed operation equations can be written and used to define the dynamic response character- 
istics of the controlled system. It also is shown how these linear mixed operation equations lead, in an 
automatable manner, directly to a set of finite difference equations which are in a format compatible 
with follow-on time and frequency domain analysis methods. 
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TIME AND FREQUENCY DOMAIN ANALYSIS OF SAMPLED DATA CONTROLLERS 

VIA MIXED OPERATION EQUATIONS 


Harold P. Frisch 

Goddard Space Flight Center 
Greenbelt, Maryland 


INTRODUCTION 

The large space systems technology (LSST) program of the National Aeronautics and Space Adminis- 
tration (NASA) is attempting to develop an interactive integrated analysis capability (IAC). Its objec- 
tive is to provide engineers with an interdisciplinary (structures-thermal-controls) analysis capability 
supportive of static, time and frequency domain design and performance evaluation methods. To do 
this, the IAC program intends to provide a time-wise efficient and user-friendly means for accessing a 
sequence of analysis tools and for respectively obtaining the necessary input data via a series of estab- 
lished interdisciplinary data flow paths. The immediate goal is to provide a coupled structures-thermal- 
controls analysis capability compatible with the inherent assumptions and limitations associated with 
current widely used state-of-the-art general purpose analysis programs such as NASTRAN, DISCOS, 
SINDA and TRASYS.* 

This stated goal implies that data must be derived and stored in computer memory in such a manner 
that it can flow between analysis modules which view intrinsically identical items of data differently. 
Today, the prime impediment to this interdisciplinary flow of data is the infinite number of ways that 
requested data can be generated and provided to the requestor. In a working environment, this prob- 
lem is usually compounded by the inability of data requestors to cross interdisciplinary semantics 
barriers and to define precisely what is being requested. The problem is particularly acute when data 
must be obtained from sources outside of one’s immediate working environment. The net effect of 
these data flow barriers is a drop in analyst productivity while assumptions and limitations associated 
with data delivered in an unfamiliar data format are investigated. 

One technique which can be used to overcome this real world problem is to define neutral formats for 
data. Conceptually, once data is in a neutral format it can, with the aid of a processor, flow anywhere. 
Since both the neutral format and the desired format are known a priori, the development of the pro- 
cessor which creates the linkage becomes a routine task. In order to minimize the number of proces- 
sors required, the choice of a neutral format converges on the answer to the question: What is the 
most commonly used, readily obtainable format for requested data? For example, in the structural 
analysis discipline NASTRAN is probably the most widely used general purpose analysis program. 
Its bulk data deck would be a good choice as the neutral format for structures data. 

'I'he choice of a neutral format for thermal analysis is complicated by the fact that there are two 
competing methods of analysis. Thermal engineers having a strong structures background usually lean 
toward the finite element approach and use either the NASTRAN or SPAR thermal analyzers. Thermal 


*See Program Acronym/ Availability List, page 18. 
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engineers having limited or no structures background usually lean toward the finite difference ap- 
proach and use the programs SINDA and TRASYS. Consequently, a bilevel neutral format will most 
probably be required for this discipline. Level one would contain the finite element neutral format 
while level two would contain the finite difference neutral format. If NASTRAN and SINDA bulk 
data decks are selected as neutral formats, then it is possible to generate an automated procedure 
which will read level 1 NASTRAN bulk data and generate from it level 2 SINDA bulk data. The reverse 
probably never will be possible. 

The choice of a single or multilevel neutral format for controls data is not at all straightforward. If the 
assumption could be made, which it cannot, that controllers are continuous system devices, then the 
input data for any of the continuous system modeling programs such as ASCL, CSMP, CSSL, DARE, 
MODEL, etc., would be candidates for a neutral format for controls data. Since most present-day 
spacecraft controllers are sampled data rather than continuous devices, new modeling programs must 
be developed or old ones extended. This is required to adequately account for, in a simultaneous 
manner, the continuous time response characteristics of the plant and the digital response character- 
istics of controller components. 

The objective of this paper is to present a framework for a neutral format which can be used for spec- 
ification of a full range of sampled data controllers. It is conceptually simple in that it parallels the 
block diagram approach conventionally used for control system documentation. Furthermore it can 
be computationally automated in such a manner that a clear distinction is made between state vari- 
ables which may be labeled as continuous time or analog signals and those which may be labeled as 
digital signals. This clear distinction allows for the development of integration algorithms which utilize 
the discontinuous nature of digital signal state variables to improve run-time performance. When the 
system is linear, this framework also leads directly to linear constant coefficient finite difference 
equations compatible with stability and frequency domain analysis methods. 

MIXED OPERATION EQUATIONS 

The concepts presented here stem largely from ideas derived from the work of Kalman and Bertram, 
1959 (Reference 4), Zadeh and Desoer, 1963 (Reference 11), and National Aeronautics and Space 
Administration/Goddard Space Flight Center (NASA/GSFC) Solar Maximum Mission (SMM) pro- 
ject support activities, Donohue et al, 1978 (Reference 2). 

Consider the problem of developing a simulation model for a linear continuous plant controlled by a 
linear sampled data control system. If this is done on an analog computer with digital components, 
then the plant equations are implemented by a set of summers, amplifiers and integrators, while the 
controller equations are implemented by a set of zero order holds, digital delays and some more 
summers, amplifiers and integrators. If this is to be done on a digital computer, the problem is not 
as straightforward due to the fact that zero order holds and digital delayors cannot be characterized by 
linear ordinary differential equations. Furthermore, if one blindly proceeds with conventional numer- 
ical integration methods, the net outcome is frequently a computationally costly solution. 

To arrive at a mathematical formulation compatible with digital simulation methods, we first look 
through the eyes of a mathematician at the bag of hardware required for the analog simulation. 
Summers and amplifiers are memoryless objects, while integrators, zero order holds and delayors are 
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finite memory objects. Carrying this a step further, we can view the integrator as a piece of hard- 
ware which implements the mathematical rules of integration; similarly, it is just as valid to view the 
piece of hardware called a zero order hold as the implementation of the mathematical rules of zero 
order hold and to view the piece of hardware called a digital delayor as the implementation of the 
mathematical rules of digital delay. 

At this point, drop the notion of hardware and recognize that new mathematical operations have 
been defined, i.e., the mathematical operations of zero order hold and that of digital delay. These 
readily lend themselves to block diagram representation and can be pictorially represented as 



where (•) dot calls for the mathematical operation of integration, (□) box calls for zero order hold 
and (A) triangle calls for digital delay. 

If several zero order holds having different sampling rates and digital delayors having different delay 
times are required, the notation immediately is expandable by introducing new symbols such as 
(O) circle, (O) hexagon, etc., or more simply by a numbering sequence within the symbols. 

How this ties in with an ability to define sampled data controllers in terms of mixed operation equa- 
tions may not yet be completely obvious. Consider, therefore, the following problem: Define the 
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mixed operation equations for the single axis control of a unit mass via sampled position and com- 
puted rate feedback through a computational delay not equal to the sampling period. The block 
diagram representation follows: 


r f 



s 


ZO H 


*3 



In this figure the following notation is used; X ; = plant rate, X 2 = plant position, X 3 = sampled 
position, (X 3 -X 4 )/At = computed rate, At = sampling period, K = position gain, D = rate gain, X s = 
control force, R F = external force input, R s = external position input. It immediately follows that 
the linear mixed operation equations which define this problem are: 




( 1 ) 
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The important point to note in this example is the ease with which the linear constant coefficient 
mixed operation equations have been defined directly from the block diagram. The constant coeffi- 
cient matrix, along with the rules associated with the mathematical operations used, completely de- 
fines the problem and can be readily implemented in a format compatible with digital computation. 

The objective now is to define the steps required to go from linear mixed operation equations to linear 
finite difference equations. Once these are in hand, follow-on frequency domain studies can be car- 
ried out. If linear time domain response is required, the finite difference equations can be coded and 
time response obtained without recourse to conventional numerical integration methods. If nonlinear 
time domain response is required, a framework for developing special purpose integration algorithms 
which account for the digital signal characteristics of certainstate variables is provided. In application, 
these numerical algorithms must often be developed to achieve timewise efficient simulation pro- 
grams. 

In the next section, the concept of mixed operation equations is presented in a more generalized 
framework and certain limitations imposed by the demand that the end product be computationally 
efficient are discussed. 

GENERALIZED FORMULATION, ASSUMPTIONS AND LIMITATIONS 

Kalman and Bertram (Reference 4) define several different types of sampling systems, e.g., synchro- 
nous, nonsynchronous, multirate, random, etc. There is no intention here to cover all possible cases. 
This formulation assumes that a single clock, contained within the system, controls all sampling. 
The clock is assumed to provide a never-ending sequence of synchronous clock pulses, the period 
between each clock pulse being exactly At. All sampling is carried out by the hardware items hereto- 
fore referred to as zero order holds. Each zero order hold is assumed to be synchronous and to have 
its sampling period defined as an integer multiple of the clock pulse period At. Zero order holds in 
any simulation may have different sampling periods; however, if frequency domain analysis is the 
goal, there must be a periodic pattern to the sampling. If not, frequency domain analysis has no 
meaning. 

In application, the penalty on generality paid for by restricting the sampling periods to be integer 
multiples of At is minimal. This is adequate for most linear and quasi-linear system design and per- 
formance evaluation studies, and it leads to mathematical requirements that can be implemented so 
as to yield timewise efficient simulation models. 

The need to incorporate the effects of one or more delays in a reasonably general format presents 
problems, if one demands that the end product be useful for general application. Delays can be placed 
into three general categories: 1) digital delay with time constant equal to an integer multiple of 
At, 2) digital delay with time constant not equal to an integer multiple of At, and 3) analog or 
transport delays. 

Analog delays are not considered here for the simple reason that the author has not been able to 
find a method which he deems to be computationally practical for timewise efficient simulations. 
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Digital delays are considered and are implemented by the use of what will be referred to here as unit 
delayors and fractional delayors subject to the following assumptions and limitations: 1) a unit 
delayor delays a digital signal, which is constant during the clock pulse period At, by the amount 
At, and 2) a fractional delayor delays a digital signal, which is constant during the clock pulse period 
At, by the amount 7At where 0 < y < 1 . These limitations are incorporated to prohibit the modeling 
of a delayed analog signal by a series connection of several fractional delayors, and to set up a frame- 
work which will lead to tractable computational requirements. 

Consider a general multivariable sampled data control system subject to the above assumptions and 
limitations, and let: 

W — column matrix containing state variables associated with all integrators 

X — column matrix containing state variables associated with all zero order holds 

Y — column matrix containing state variables associated with all unit delayors 

Z — column matrix containing state variables associated with all fractional delayors 

• — symbol for integration 

□ - symbol for zero order hold; subpartitioned into CD , 0 , . . . etc., if many different zero 
order holds are used 

A — symbol for unit delay 

O - symbol for fractional delay, subpartitioned into <D , <D , . . . etc., if many different 
fractional delays are used 

Based upon the above comments and notation, the general form for a system of linear mixed opera- 
tion equations is 
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where 


A,B — partitions with nonzero elements 
0 — partitions with all zero elements 

R - column matrix of reference inputs assumed to be continuous functions of time; these 
cannot feed directly into a delayor. 

It is also possible to define systems of nonlinear mixed operation equations. These, in general, will 
not painlessly lead to finite difference equations. However, at times they do lead to special-purpose 
integration algorithms which are superior to those specifically developed for systems of first order 
nonlinear ordinary differential equations. The user is cautioned that relative superiority must be 
judged on a problem-by-problem basis and is strongly dependent upon the talents of the analyst/ 
programmer. 

MIXED MATHEMATICAL OPERATIONS 

The system of linear mixed operation equations as presented is relatively easy to obtain directly 
from a block diagram definition of a controller. However, in order for these to be useful for follow- 
on analysis, they must lead to a set of either linear ordinary differential or linear finite difference 
equations in an automatable manner. The objective of this section is to outline the methodology to 
be used to automate the procedure by which system state at clock time (k+l)At can be expressed 
as a linear function of system state and input at clock time kAt, k = 0,1,2,.... This is the desired 
set of first order linear finite difference equations which, via use of the z-transform, will lead directly 
into follow-on frequency domain analysis methods. 

Let W(k), X(k), Y(k), Z(k) define system state and let R(k) define system input at clock time kAt. 
To define system state at time (k+l)At, it is necessary to know the time dependency of all system 
state variables and inputs during the At time interval 


kAt < t < (k+1) At 


This is known via the modeling restrictions imposed, and from the mathematical rules associated 
with the mixed operation equations given below as: 

• Integration: The matrix equation which defines all integrator state variables takes on the 
general form 


W(t) = Aj j W(t) + U(t) 


(3) 
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where U(t) is everything in row 1 of equation 2 not explicitly stated herein. Let W(k) define a full 
set of initial conditions for equation 3. It follows from the mathematical rules associated with the 
operation of integration that the state of all integrator state variables at any time t in the time interval 


kAt < t < (k+l)At 


is given by 


W(t) 


= e^i ' ^^W(k) + f e^i i ^ U(r) d r 

KAt 


(4) 


where it is immediately pointed out that numerical algorithms for evaluating the matrix exponential, 
and an integral involving the matrix exponential, will be required. From the standpoint of computa- 
tional mathematics, the numerics associated with the evaluation of equation 4 will be practical only 
if the components of U(r) are digital and delayed digital signals with no more than one jump dis- 
continuity in the integration interval. This is the reason for the modeling restrictions which prohibit 
analog delayors and series connections of fractional delayors. 

• Zero Order Hold: The matrix equation which defines all zero order hold state variables takes on 
the general form 


□ 

X(t) = U(t) 


(5) 


where U(t) now is everything on the right-hand side of equation 2 in row 2. A zero order hold 
samples the incoming signal at the sampling instant, it updates its output signal an instant later 
and then holds that constant value until it is subsequently updated an instant after the next sampling 
instant. The output signal of a zero order hold is frequently referred to herein as a digital signal. 
The above descriptive definition of what a zero order hold does is encapsulated in the mathematical 
solution to the zero order hold equation (5) given by 


X(t) = U(k) 


for 


kAt < t < (k+l)At 


( 6 ) 



If zero order holds with different sampling rates must be studied, one cannot blindly apply equation 6. 
For example, assume a zero order hold with sample period nAt. A trivial notation extension allows 
one to symbolically express its equation as 


X(t) = U(t) (7) 


If clock time kAt is a sampling instant for this zero order hold, then the solution to equation 7 
becomes 


for 


X(t)=U(k) 
kAt < t < (k+n)At 


( 8 ) 


• Unit Delay: The matrix equation which defines all unit delay state variables takes on the 
general form 


A 

Y(t) = U(t) (9) 


where U(t) now is everything on the right-hand side of equation 2 in row 3. A unit delay reads the 
incoming digital signal at the sampling instant kAt, and then outputs the state it was at, at kAt, until 
an instant before the next sampling instant at which time it updates its output signal. Mathema- 
tically, this is encapsulated in the solution to the unit delay equation which is given by 


Y(t) = 


Y(k) 

U(k) 


for kAt < t < (k+ 1 ) At 
for t=(k+l)At 


( 10 ) 


® Fractional Delay: The matrix equation which defines all fractional delay state variables takes 
on the general form 


O 

Z(t) = U(t) (11) 


where U(t) now is everything on the right-hand side of equation 2 in row 4. A fractional delay reads 
the incoming digital signal at the sampling instant kAt. It then outputs the value of its state, at kAt, 
until the time (k+ 7 )At at which time it updates its output signal. Its output signal is referred to 
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herein as a delayed digital signal. Mathematically, this is encapsulated in the solution to the fractional 
delay equation, i.e., 


Z(t) = 


' Z(k) 

< 

' U(k) 


kAt < t < (k+7)At 
(k+y)At<t<(k+l)At 


( 12 ) 


where 0 < 7 < 1. If fractional delayors with different delay times are required for simulation, sub- 
scripts would be required on the 7 ’s and in the symbol O in a manner similar to that used in the 
zero order hold discussion. 

SOLUTION OF MIXED OPERATION EQUATIONS 

The most general form for linear mixed operation equations considered in the context of this paper 
is presented in equation 2, This can be expressed as a set of first order finite difference equations by 
direct application of the closed-form solutions for each operation considered in the previous para- 
graphs. In keeping with the author’s intent to communicate to the reader, the solution presented 
herein assumes that all zero order holds have sampling period At and all fractional delays have delay 
time 7 At. The solution for the multizero order hold, multifractional delay situation requires addi- 
tional notation which on the bottom line is more confusing than illuminating; it will therefore not 
be presented here. 

The row 1 equation of equation 2 defines the integrator state variables and is given by 

W(t) = A, , W(t) + A, 2 X(t) + Aj 3 Y(t) + A 14 Z(t) + B, R(t) (13) 

In order to concisely define the solution to this equation let 


A .At 

F(At) = e 11 


At 


Aj , (At-r) 


H(At) = f e “ dr 


At A (At-r) 

H [ ( 1 - 7 ) At ] = / e 11 dr 


/ /v 


7At 


(14) 


(15) 


(16) 
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and note that 


H(7At) = H(At)-H[(l-7)At] 


(17) 


Direct application of equation 4 evaluated at t = (k+l)At along with equations 6, 8, 10 and 12 
leads to the finite difference equation 


W(k+1) = F(At) W(k) + A, 2 H(At) X(k+1) 

+A H(At) Y(k) 

13 (18) 
+ A i4 |H( 7 At)Z(k) + H[(l- 7 )At] Z(k+l)| 

+Bj H(At) R(k) 


where the assumption is made that over the At time period the input R(t) is constant and equal to 
R(k). Furthermore, note the use of X(k+1) and Z(k+1) in the right-hand side of equation 18. This 
is a statement of the fact that over the interval of integration, these state variables are defined by 
their state at (k+l)At. 

The row 2 equation of equation 2 defines the zero order hold state variables. Direct application 
of equation 6 leads to the finite difference equation 


X(k+ 1 ) = A 2 , W(k) + A 2 2 X(k) + A 2 3 Y(k) + A 2 4 Z(k) + B 2 R(k) ( 1 9) 


The row 3 equation of equation 2 defines the unit delay state variables. Direct application of 
equation 1 0 leads to the finite difference equation 


Y(k+1) = A 32 X(k) + A 33 Y(k) +A 34 Z(k) (20) 
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The row 4 equation of equation 2 defines the fractional delay state variables. Direct application of 
equation 12 leads to the finite difference equation 


Z(k+1) = A 42 X(k+l) + A 43 Y(k+l) (21) 

Equations 18, 19, 20 and 21 can now be combined to yield the following matrix of simultaneous 
first order finite difference equations 
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Direct use of the inverse of the coefficient matrix on the left-hand side of the equation, which is 


1 (A 12 H(At) {A 14 H[(l- 7 )At]A 43 | { A, 4 H[(l- 7 )At] \ 

+A 14 H[(l-r)At]A 42 } 


0 

0 


4 2 


4 3 


(23) 


0 

1 


yields the final set of finite difference equations which are in a format compatible with follow- 
on analysis methods, i.e., 


( W(k+1) ^ 

| 

~ D m 

D 12 

D 13 

D lf 

) X(k+1) | 

s. = 

D 2 1 

D 2 2 

D 2 3 

D 2 4 

I Y(k+1) 

r 

° 31 

D 3 2 

D 3 3 

d 34 

(z(k+l) ^ 


_ D <‘ 

D 4 2 

D 4 3 

D 4 4 




DIGITAL COMPUTATIONAL REQUIREMENTS 

Wilkinson (Reference 10) reminds us that “attractive mathematics does not protect one from the 
rigors of digital computation.” It is therefore important to determine whether or not the foregoing 
presentation has led to a computationally practical methodology for the particular case worked 
out and for the multizero order hold, multifractional delay case mentioned. 

The following sequence of steps is required to go from problem definition to the desired set of finite 
difference equations to be used for either time or frequency domain follow-on analysis. 
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Step 1. Set up the linear mixed operation equations in the partitioned matrix format of equation 
2. This is straightforward and can be done by hand or problem-unique computer code. 

Step 2. Determine the elements in the coefficient matrices on the right- and left-hand side of 
equation 22. The major numerical problem here is in the evaluation of the matrix exponential 
F(At) and the integrals involving the matrix exponential, namely, H(At) and H(7At). These matrices 
can be efficiently and accurately evaluated by use of subroutine PADE developed by Van Loan and 
discussed in Reference 8, discussed and published in Reference 9. 

Step 3. Determine the elements in the coefficient matrices of equation 24. These are obtainable 
by a simple matrix multiply, since the required matrix inverse is definable in closed form. 

FREQUENCY DOMAIN ANALYSIS 

If frequency domain analysis is the objective, a periodic pattern to the sampling must exist. If there 
exists a periodic pattern to the sampling, then it is possible via successive matrix multiplications to 
arrive at an equation of the form 


W(k+N) 

X(k+N) 

Y(k+N) 

Z(k+N) 


* 1 , 

*12 

*13 

*14 

/ W(k)\ 

v i 

*2 1 

*2 2 

*2 3 

*24 

) X(k) 1 

V 2 





\ / + 


*3 1 

*32 

*3 3 

*34 

) Y(k) l 


*4 1 

*4 2 

*4 3 

*4 4 

yz(k) J 

\ 



(25) 


where NAt is the sampling pattern period and the elements of the coefficient matrix are constant 
and independent of both k and N. Kalman and Bertram (Reference 4) discuss this and refer to this 
matrix as a stationary transition matrix. It should be noted that in equation 25 the standard stability 
analysis assumption that input R(k) is constant over the pattern period NAt has been made; if not 
made, the resultant stationary finite difference equation is unmanageable for follow-on frequency 
domain analysis. 
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In order to transform time domain equation 25 into a frequency domain equation of standard form, 
let AT = nAt define the sampling pattern period and define the z-transform of the arbitrary discrete 
time function f(nAT) n=o, 1, 2, ...as 


f(z) = ^2 f ( nAT > z ' n (26) 

n=0 


Recognize that W(nAT), X(nAT),... and R(nAT) are discrete time functions, utilize equations 26 to 
transform equation 25 and obtain the desired set of standard form frequency domain equations, i.e., 


[zU- )//] 



= [v] R(z) 


(27) 


where "fl is the unit diagonal matrix and the matrices \p and v are the coefficient matrices defined in 
equation 25. Procedural details required for the above transformation can be found in Reference 3 
(Jury, 1964). 

Tire reader should note the author’s intentional avoidance of the mention of the Laplace transform 
variable s in the z-transform definition used. From the standpoint of computational mathematics, 
it is the author’s opinion that numerical computation should be carried out with the z-domain equa- 
tions provided in equation 27. Once z-domain results are computationally obtained, they can always 
be mapped into: 

The s-domain by the mapping equation 


z = e sAt 


The w-domain by the mapping equation 


z = 


1 + w 
1 - w 
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or the £,-domain by the mapping equation 


C = 


z - 1 


At 


This final mapping step should be viewed as a necessity since results must be presented in a format 
compatible with the desires of control design engineers. This computational approach avoids truncated 
series approximations, and the resultant numerical computation pitfalls which can be encountered, 
when computation is done with s-domain equations. 

Determination of closed loop system stability can be made directly from the eigenvalues of the 
coefficient matrix in equation 25; these may be computed via application of the EISPACIC library 
of matrix eigensystem routines published in Reference 7 (Smith, et al., 1976). If all eigenvalues are 
contained within the unit circle in the complex z-domain, then the system is stable. If any eigenvalue 
is outside of the unit circle, the system is unstable. 

Determination of margins of stability is usually accomplished with opened, partially opened and 
closed loop z-domain transfer functions. A methodology for generating any of these type transfer 
functions which is based upon a state variable formulation of system response is presented in Refer- 
ence 1 (Bodley, et al., 1978) and implemented in the DISCOS program defined therein. Even though 
the DISCOS method is specifically for continuous time systems, and hence s-domain transfer func- 
tions, slight modifications can be introduced which make it applicable for z-domain transfer functions. 

The DISCOS method opens control feedback loops by selectively setting elements to zero in the 
coefficient matrix of the set of first order ordinary differential equations used to define system 
response. This resultant coefficient matrix is the coefficient matrix for an opened loop system, and 
resultant closed loop transfer functions are in actuality the sought-after open loop transfer functions 
of the original unabridged system. To apply the method to a sampled data system, it is important to 
recognize that feedback loops must be opened by setting to zero elements in the coefficient matrix 
of the mixed operation equations. One then proceeds directly as defined herein to obtain the coeffi- 
cient matrix for the finite difference equations, which define open loop system response, and then on 
to z-domain transfer function analysis. 

TIME DOMAIN ANALYSIS 

If time domain analysis is the objective, the concept of mixed operation equations, at times, can lead 
to tailor-made integration methods which, in particular applications, can reduce computation run 
time by orders of magnitude. 

For example, if the matrix in equation 2 is constant and input R is constant over each At 
sampling interval, classic numerical integration methods for the differential equations can be com- 
pletely avoided. One needs only pass through the matrix exponentiation routine PADE once, save 
the results, and then compute time domain response directly from the finite difference equations 
provided by equation 24. If the matrix A ( is slowly varying relative to At, it often will be found 
to be more efficient to call up the matrix exponentiation routine periodically to update data for 
equation 24, rather than blindly use classical numerical integration methods. 
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The reader should recognize that numerical integration methods for differential equations are the 
subject of active research; however, most of the current work being done is concentrated on what in 
the literature is referred to as stiff ordinary differential equations (Shampine, et ah, 1976 and 1979, 
References 5 and 6, and Gear, 1981, Reference 12). Methods for integrating equations of the type 
encountered in the sampled data control of oscillatory plants, as defined herein, is an area of research 
not being addressed in the open literature today. Hopefully, the ideas presented in this paper will 
stimulate research into the field of numerical integration of mixed operation equations. This is an 
important class of equations for which timewise efficient methods for numerical solution are now 
developed on a problem-by-problem basis. 

CONCLUSION 

The methodology presented here is a generalization of methods successfully used at NASA/GSFC 
during a control system stability and performance analysis of several modes of operation for the 
SMM spacecraft. In that application equation 24 was transformed by use of the z-transform, so that 
stability margins could be determined in the frequency domain. As a benchmark; in that application, 
20 different modes of spacecraft operation were analyzed, each requiring the regeneration of the 
coefficient matrix in equation 24. Seventeen state variables were required to define the plant and 
controller, and all four mathematical operations discussed herein were utilized. Upper and lower 
gain margins were obtained in a straightforward iterative manner by incrementing feedback gain 
and checking for the gain values which put system roots (eigenvalues) on the threshold of instability, 
i.e., out of the z-domain unit circle. Once debugged, the entire job was run in 1.6 minutes CPU on the 
IBM 360-91. Results were compared with those obtained by Donohue, et al., 1978, (Reference 2), 
using more conventional s-domain methods. Agreement in all cases was near perfect. 

Donohue, Hager, McGlew and Zimmerman also carried out time domain performance studies for the 
SMM spacecraft. The tailor-made integration method used in their study contributed significantly 
to the ideas used and presented here in a generalized context. Their experience conclusively demon- 
strated that for certain real application problems, dramatic (orders of magnitude) improvements in 
computer run-time could be obtained if one avoided the blind application of conventional numerical 
integration methods for sampled data control equations. 
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